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Several  multigrid  schemes  are  considered  for  the  numerical  computation  of  viscous 
hypersonic  flows.  For  each  scheme,  the  basic  solution  algorithm  employs  upwind  spatial 
discretization  with  explicit  multistage  time  stepping.  Two-level  versions  of  the  various 
multigrid  algorithms  are  applied  to  the  two-dimensional  advection  equation,  and  Fourier 
analysis  is  used  to  determine  their  damping  properties.  The  capabilities  of  the  multigrid 
methods  are  assessed  by  solving  three  different  hypersonic  flow  problems.  Some  new 
multigrid  schemes  based  on  semicoarsening  strategies  are  shown  to  be  quite  effective 
in  relieving  the  stiffness  caused  by  the  high-aspect-ratio  cells  required  to  resolve  high 
Reynolds  number  flows.  These  schemes  exhibit  good  convergence  rates  foi  Reynolds 
numbers  up  to  200  x  10*^  and  Mach  numbers  up  to  25. 
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NASA  Co  .tract  No.  NASl-18605  while  the  first  author  was  in  residence  at  the  Institute  for  Computer 
Applications  in  Science  and  Engineering  (ICASE),  NASA  Langley  Research  Center,  Hampton,  VA 
23665. 
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1  Introduction 


Over  the  past  few  years  the  need  for  efficient  numerical  methods  to  solve  the  equations 
governing  hypersonic  viscous  flows  has  become  very  obvious.  Mostly,  flow  solvers  used  in 
current  aerospace  programs,  such  as  X-30  or  SAENGER,  exhibit  slow  convergence  towards 
the  desired  steady-state  solutions,  which  leads  to  high  computer  costs,  long  turn-around 
times,  and  a  slowdown  in  the  efforts  to  design  these  vehicles.  The  reason  for  this  is  the 
appearance  of  flow  phenomena  with  very  different  scales  and  with  highly  nonlinear  behav¬ 
ior.  We  mention  here  the  laminar  and  turbulent  boundary  layers  at  very  high  Reynolds 
numbers  and  their  interactions  with  shocks  and  slip  lines,  and  furthermore,  shock/shock 
interactions  which  generate  complex  flow  fields.  Many  numerical  techniques  which  were 
developed  to  assist  convergence  of  subsonic  or  traiisonic  flow  calculations  are  found  to  be 
inappropriate  for  hypersonic  flow  applications.  For  example,  the  time  step  of  many  explicit 
and  implicit  schemes,  which  allows  the  transient  behavior  of  strongly  nonlinear  hypersonic 
flow  phenomena  to  be  captured,  is  highly  limited.  Consequently,  thousands  of  time  steps 
are  needed  to  converge  the  thin  boundary  layers. 

A  particular  method  which  was  successfully  developed  to  accelerate  convergence  for 
a  broad  range  of  flow  problems  at  both  subsonic  and  transonic  speeds  is  the  multigrid 
approach.  This  method,  which  uses  a  sequence  of  successively  coarser  meshes  in  ord..r  to 
propagate  disturbances  throughout  the  flow  field,  combines  nicely  with  explicit  multistage 
time-stepping  schemes  [1],  Good  convergence  rates  were  obtained  for  inviscid  flow  and, 
later  on,  for  viscous  flows  also  [2-4].  Initial  attempts  to  apply  this  promising  method  to 
hypersonic  flows  failed  fr  several  reasons.  Primarily,  the  shock  capturing  capabilities  of 
the  central- difference  scheme  used  in  [1-4]  were  found  insufficient  to  resolve  strong  shocks. 
Subsequently,  the  shock  detection  mechanism  built  into  the  central- difference  scheme  was 
improved  in  [5,6].  In  order  to  have  strong  shocks  and  slip  lines  resolved  with  fewer  com¬ 
putational  points,  the  central-difference  scheme  was  replaced  with  an  upwind-type  scheme 
in  [7-9].  Since  the  high-frequency  damping  properties  of  upwind  schemes  are  generally  less 
controllable  compared  with  central-difference  schemes,  a  variant  of  the  standard  multigrid 
approach  was  also  used  in  [8,9].  With  this  variant,  additional  coarse  meshes  are  generated 
by  semicoarsening  in  the  different  coordinate  directions.  This  strategy  was  felt  neccessary 
to  alleviate  convergence  problems  associated  with  the  high-aspect-ratio  cells  of  the  com¬ 
putational  mesh.,  An  additional  problem  encountered  is  that  very  high  temperatures  can 
occur  in  the  stagnation  region  and  near  the  surface  at  high  Mach  numbers;  and  hence, 
the  time  step  of  explicit  schemes  may  be  severely  restricted  by  the  viscous  stability  limit.; 
It  was  found  [6,8,9]  that  the  viscous  time-step  limit  can  be  removed  by  implicit  residual 
averaging. 

It  is  worthwhile  to  notice  additional  published  work  on  nuiltigrid  schemes  for  hy¬ 
personic  flows.  Decker  and  Turkel  [10]  analyzed  the  effect  of  boundary  conditions  and 
Runge-Kutta  coefficients  on  multigrid  convergence  for  hypersonic  inviscid  flows.,  Leclercq 
[11]  analyzed  two- level  multigrid  cycles  with  multistage  schemes  and  upwind  differencing 
in  one  dimension,  and  she  presented  two-dimensional  computations  of  hypersonic  flows  on 
unstructured  meshes..  As  a  means  to  remove  stiffness  associated  with  high-aspect-ratio 
cells,  Blazek  et  al.  [12]  introduced  an  upwind-biased  form  of  the  residual  smoothing  by 
which  higher  Courant  numbers  could  be  obtained.,  Thomas  [13]  used  multigrid  in  combina- 
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tion  with  third-order  upwind  difFerencing  and  implicit  approximate  factorization  schemes. 
Koren  and  Hemker  [14]  solved  the  steady  Euler  equations  with  multigrid  and  point  relax¬ 
ation  applied  as  the  smoother.  Using  damped  restriction  and  upwind  prolongation,  they 
reported  impressive  convergence  rates  for  high-speed  flows  around  blunt  bodies. 

The  present  paper  describes  recent  efforts  to  understand  and  to  improve  the  use 
of  multigrid  schemes  for  the  computation  of  hypersonic  flows.  First,  various  two- level 
multigrid  schemes  with  and  without  semicoarsening  are  introduced.  Then  we  use  Fourier 
analysis  of  the  schemes,  when  applied  to  the  two-dimensional  convection  equation,  in  order 
to  study  the  behavior  of  their  components.  For  each  multigrid  approach,  the  solver  uses 
upwind  discretization  combined  with  an  explicit  multistage  scheme.  We  next  consider  the 
numerical  solution  of  the  Navier-Stokes  equations  for  hypersonic  flows.  In  Section  5,  the 
basic  elements  of  the  flow  solver  for  these  equations  are  described.  Some  details  concerning 
the  application  of  the  time-stepping  scheme  to  fine  and  coarse  grid  problems  are  presented 
in  the  first  part  of  Section  6.  The  extension  of  the  two-level  schemes  to  multilevel  ones 
is  then  discussed.  Elements  of  multigrid  that  are  of  particular  importance  for  high-speed 
flow  computations  are  given.  In  the  results  section,  we  consider  three  different  hypersonic 
flow  problems  to  assess  the  capabilities  of  the  multigrid  schemes.  The  effect  of  stiffness, 
arising  from  coordinate  grids  with  high-aspect-ratio  cells  and  from  flow  alignment,  on 
the  performance  of  the  multigrid  methods  is  examined..  The  benefits  of  semicoarsening 
are  clearly  demonstrated.  Moreover,  with  the  semicoarsening  strategies  being  considered, 
good  convergence  rates  are  obtained  for  Reynolds  numbers  up  to  200  x  10®  and  Mach 
numbers  up  to  25. 
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2  Multigrid  Strategies 


To  set  the  stage  for  the  discussion  relating  to  multigrid  in  subsequent  sections  of  this 
paper,  we  first  briefly  describe  the  multigrid  method  and  different  execution  strategies 
that  will  be  considered.  The  multigrid  approach  is  based  on  the  Full  Approximation 
Scheme  of  Brandt  [15].  The  grid  transfer  operators  are  those  considered  by  Jameson 
[Ij.  Coarser  meshes  are  obtained  by  eliminating  alternate  mesh  points  in  each  coordinate 
direction.  Both  the  solution  and  the  residuals  are  restricted  from  fine  to  coarse  meshes. 
A  forcing  function  is  constructed  so  that  the  solution  on  a  coarse  mesh  is  driven  by  the 
residuals  collected  on  the  next  finer  mesh.  The  corrections  obtained  on  the  coarse  mesh 
are  interpolated  back  to  the  fine  mesh.  The  multigrid  schemes  investigated  within  the 
present  work  are  displayed  in  Figure  1.  Figure  1  a)  shows  a  two- level  scheme  with  full 
coarsening.  Restriction  of  the  solution  from  the  fine  mesh,  (m,n),  to  the  coarse  mesh, 
(m/2,n/2),  is  done  by  injection,  whereas  full  weighting  is  used  for  the  restriction  of  the 
residuals.  Prolongation  of  the  corrections  is  done  by  bilinear  interpolation.  Figure  1  b) 
shows  a  scheme  with  semicoarsening  in  the  different  coordinate  directions.  Again,  injection 
and  full  weighting  are  used  in  the  restriction  process.  The  corrections  obtained  on  the 
coarse  meshes  are  averaged  before  adding  them  to  the  fine  mesh.  This  is  indicated  by  the 
numbers  at  the  “up”  arrows..  Due  to  this  averaging,  half  of  the  individual  corrections  on  the 
coarse  meshes  is  lost.  It  is,  therefore,  anticipated  that  the  scheme  in  Figure  1  a^i  should  be 
computationally  more  efficient,  provided  there  is  enough  high-frequency  damping  obtained 
with  the  smoothing  scheme  of  the  fine  mesh.  In  order  to  overcome  this  deficiency  of  the 
semicoarsening  scheme,  two  more  variants  are  considered.  For  the  scheme  of  Figure  1 
c),  the  solutions  on  the  coarse  meshes  are  computed  sequentially.  Hence,  the  corrections 
obtained  on  the  (m/2,n)  mesh  can  be  used  to  update  the  (m,n/2)  mesh  before  time  stepping 
(as  indicated  by  the  horizontal  arrow).  The  sequential  update  of  the  second  coarse  mesh 
allows  the  full  amount  of  corrections  to  be  passed  up  to  the  fine  mesh.  Note  that  this 
multigrid  variant  is  not  compatible  with  the  idea  of  parallel  computations.  An  interesting 
compromise  between  schemes  of  Figures  1  b)  and  1  c)  was  suggested  by  Van  Rosendale 
[16]  (Figure  1  d)).  Here,  only  the  corrections  common  to  both  of  the  coarse  meshes, 
(m/2,n)  and  (m,n/2),  are  averaged,  whereas  the  corrections  to  the  modes  living  either  on 
(m/2,n)  or  on  (m,n/2)  are  passed  to  the  fine  mesh  in  full.  This  scheme  does  allow  parallel 
computations  for  the  coarse  meshes. 

3  Fourier  Analysis  of  the  Scalar  Advection  Equation 

A  crucial  factor  in  constructing  an  effective  multigrid  method  is  the  selection  of  a 
smoothing  or  driving  scheme.;  Local  mode  (Fourier)  analysis  is  generally  applied  to  evaluate 
possible  smoothers  on  the  basis  of  stability  and  high-frequency  damping  properties..  The 
screening  of  schemes  is  often  performed  with  a  single-grid  analysis.;  Since  a  stable  single¬ 
grid  scheme  may  not  be  stable  for  the  multigrid  process,  the  behavior  of  a  smoother  with 
a  particular  multigrid  strategy  is  needed..  In  addition,  the  multigrid  process  can  have 
a  substantial  impact  on  the  performance  of  the  multigrid  method.  In  fact,  as  we  will 
demonstrate  in  this  paper,  semicoarsening  can  provide  significant  improvement,  relative 
to  full  coarsening,  in  the  damping  of  the  multigrid,  especially  when  there  is  a  strong  mesh 
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anisotropy  due  to  high-aspect-ratio  cells. 

A  two-grid  or  two-level  multigrid  analysis  has  been  applied  by  Jameson  [1],  Mulder 
[17],  and  Leclercq  [11]  using  various  schemes  (i.e,,  multistage  time  stepping,  different  types 
of  relaxation)  for  solving  the  Euler  and  Navier-Stokes  equations  of  fluid  dynamics.  In 
[1]  Jameson  introduces  a  multilevel  uniform  analysis  and  considers  the  linear  advection 
equation  in  one  space  dimension.  This  approach  represents  a  departure  from  the  standard 
two-grid  analysis  given  by  Stiiben  and  Troitenberg  [18],  which  forms  the  basis  for  the 
analysis  used  in  [17]  and  [11].  With  the  multilevel  uniform  analysis,  fine-grid  and  coarse- 
grid  corrections  are  computed  at  all  points  of  the  fine  grid.  Then  a  nonlinear  filter  is 
applied  to  remove  the  coarse-grid  corrections  at  fine-grid  points  not  contained  in  the  coarse 
grid.  The  filtering  produces  additional  errors  in  the  form  of  a  carrier  wave  with  a  frequency 
depending  on  the  fine-mesh  spacing.  This  analysis  does  not  allow  for  the  coupling  (aliasing) 
effects  due  to  the  restriction  operator  (fine  to  coarse  grid  transfer  operator)  in  the  multigrid 
metiiod..  However,  it  does  offer  the  advantages  of  simplicity  and  application  to  more  than 
two-level  schemes.  Thus,  it  allows  the  rapid  comparison  of  multigrid  algorithms.  If  a 
multigrid  method  is  unstable  or  inefficient  according  to  this  analysis,  then  it  is  probably 
not  a  reasonable  scheme. 

In  this  section,  we  consider  the  scalar  two-dimensional  advection  equation.  The  mul¬ 
tilevel  uniform  analysis  of  [1]  is  extended  to  two  space  dimen3ions  and  applied  to  full 
coarsening  and  semicoarsening  strategies. 

Consider  an  initial  value  problem  governed  by  the  scalar  advection  equation 

~w{x,iy,t)  +  Cw(x,yJ)  =  0,  {x,xj)  6  (3.1) 


where  the  linear  operator  is 


£  = 


d  ,  a 


a  >  0,  6  >  0, 


the  domain  C  and  f  €  Sf"'".  Assume  a  periodic  boundary  condition  for  the  scalar 
function  w{x,y,t).  Let  1)  =  {(;r,(/) 0  <  x  <  1,  0  <  i/  <  1}.;  Define  a  fine  grid  Gf  and 
coarse  grids  Gc,i  {I  —  1,2)  tha'  cover  the  domain  D  such  that  Gc,i  C  Gf.  We  generate 
grids  Gc,t  by  eliminating  every  other  me.sh  line  of  G  f  in  one  or  both  coordinate  directions. 
First  we  describe  the  fine-grid  discrete  problem.  Let  the  grid  G /  contain  m/  xuf  cells  and 
have  uniform  spacings  Ax/  and  Ay/.  Let  the  discrete  function  (t(^/)t,j  reside  at  the  Gf 
mesh  point  (iAx/,  j  Ay/).  At  each  point  (i,  j),  we  consider  a  corresponding  cell  {Cf),j  with 
corners  at  (/-l/2,j-l/2),  («  + 1/2,/ -1/2),  (i-f  l/2,/-t-i/2),  and  (i-l/2,/  +  l/2)..  Suppose 
we  approximate  the  spatial  derivatives  of  (3.1)  with  first-order  upwind  differencing.  Then, 
we  obtain 


^tfj^UvfUj  =  -N^  [{wf)l  -  {«>/);-!, ,1  -  [(.e/)"^  -  {wf)\\^_,]  ,  (3.2) 


where  the  Courant  numbers  are 


=  A^(T/,  N,,  =  X,,af 
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=  aAy/,  A,  =  bAxj,  a/  = 

The  superscript  n  denotes  the  time  level  n At.  If  we  estimate  the  time  step  At /  by 


Atf  = 


_  NAilf 


aAyf  +  hAxf' 


it  follows  that 


aN 

a  +  &A/  ’ 


Nr, 


bNAf 
a +  bAf' 


where  N  is  the  Courant- Friedrichs- Lewy  (CFL)  number  for  an  explicit  time-stepping 
scheme,  and  A/  is  the  cell  aspect  ratio  AxffAy/.  Taking  a  Fourier  transform  of  (3.2),  we 
obtain 


where 


f(6)  =  (1  -  cosB)  -1-  isin^, 
and  1  =  a/^.  The  transformed  discrete  function  is 

m/— 1  n/— 1 

1 1=0  ji=0 

where  the  phase  angles  0^  and  dn  are  given  by 

=  27r-^,  =  27r— , 

m/  Uf 


with  wave  numbers 


Equation  (3.5)  can  be  rewritten  as 

6w,  =  (»,)"+■  -  (w,r  =  (3.7) 

with  Z  being  the  Fourier  symbol  for  the  difference  operator.  Consider  the  explicit  p-stage 
scheme 

=  (<./)" 

.'=l,2,--,p  (3.8) 


o 


where  i2/  is  a  residual  function  defined  as 


Rf  =  AQfCfWf. 

Using  this  scheme  and  (3.7),  one  can  always  represent  the  change  in  w/  by 

and  the  symbol  of  the  time-stepping  operator  F  is  given  by 

=  ^P~  OlpOip-lZ  -1-  0lp0lp-x(\p-2Z^ - 


(3.9) 

(3.10) 


(3.11) 


In  (3.11)  we  have  assumed  that  both  the  convective  and  dissipative  parts  of  the  symbol 
Z  are  evaluated  on  each  stage  of  the  time-stepping  scheme.  This  scheme  is  a  member  of  a 
general  class  of  (p,  q)  schemes..  The  p  refers  to  the  number  of  stages,  and  q  designates  the 
number  of  evaluations  of  the  dissipative  contribution..  Suppose  we  have  a  (p,  2)  scheme. 
Let 

Zr  =  Re{Z),  Zi  —  Im{Z). 

Then,  if  p  >  3,  the  F  of  (3.11)  is  replaced  by 

F  =  Op  -  Oip{a\ZR  -f  (Xp-iZi)Z'i  -I-  Oipap-\{a\ZR  -f  oip^2Zi)Z] - 

-  (-l)'-'(apa,-i  •  •  ajXo.Zfl  +  a2Z,)Zr’  (3-12) 

+  (— 1)^  1  ■  ■  ■ 

In  this  paper,  we  consider  a  (5, 3)  scheme  where  the  dissipative  evaluations  are  weighted. 
For  this  scheme,  the  symbol  of  the  residual  function  corresponding  to  the  (A;  -j-  l)th  stage 
is  written  as 


^(^■+1)  ^  ^-1 


k 

1=0 


k 

Y^'rki  =  i> 


1=0 


The  weighting  faciors  are  as  follows: 


7oo 

=  1, 

7io 

=  1, 

711 

=  0, 

720 

=  r3, 

721 

=  0, 

722 

=  73. 

730 

=  r3, 

731 

=  0, 

732 

=  73. 

733 

=  0, 

740 

=  r3F,3, 

741 

=  0, 

742 

=  7.3^5, 

743 

=  0, 

(3.13) 


(3.14) 


where  Fa  =  (1  —  73),  Fs  =  (1  —  75),  73  =  0.5C,  and  75  =  0.44.  The  symbol  of  the 
time-stepping  operator  is  given  by 


F  =  05  [1  —  04^1(1  —  oisZi)  —  a^Z:iZ\{oi^Z2Zi  —  ■j^Zr)  —  F573Z3Z/i] ,  (3.15) 
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where 

Zi  =  2/  +  75^ft, 

Z2  =  Zj+  ^^Zr, 

Zz  ~  a2{l  —  oiiZi). 

One  can  extend  the  stability  range  of  the  explicit  scheme  of  (3.8)  with  implicit  residual 
smoothing.  For  two  dimensions,  the  residual  smoothing  can  be  applied  in  the  form 

(1  -  ;9«V<Aj)(l  -  =  nf,  (3.16) 

where  the  residual  is  defined  using  (3.9)  as 

(n,f>  =  (3.17) 


A  and  V  are  the  usual  forward  and  backward  difference  operators,  and  I  =  1,  •  •  •  ,p,  with  p 
being  the  number  of  stages  in  the  scheme.  The  variable  coefficients  {3^  and  are  defined 
as  follows: 


=  max 


1 

4 


/jV  1 

\N*  l  +  il>rrf^ 


(3.18) 


=  max 


1 

4 


IL  1 

N*  1  +  ^r-' 


,0 


with  N/N*  the  ratio  of  the  Courant  number  for  the  smoothed  scheme  to  that  of  the 
unsmoothed  scheme,  and  =  A,,/A^,  A  reasonable  value  for  the  parameter  ^  is  0.125. 
In  practice,  this  implicit  procedure  allows  the  explicit  stability  limit  to  be  increased  by  a 
factor  of  2  to  3.  For  additional  details  concerning  implicit  residual  smoothing,  see  [3,5]. 

If  residual  smoothing  is  applied  on  each  stage  of  the  time-stepping  scheme,  then  the 
symbol  Z{6(,dq)  appearing  in  (3.10)  and  (3.11)  is  replaced  with 


/(«{,«,)  =  s(«<.«,)z(«5,«,),  s(«j,6,)  =  rj-'r;',  (3.19) 


where 

=  H-  2/?^(l  -  cos^4),  r,,  =  1  -f  2/9,,(l  -  cos^,,). 

Before  considering  the  discrete  problem  on  some  coarse  grid  Gc,i,  we  define  the  re¬ 
striction  operators  and  their  corresponding  Fourier  symbols.  Assume  that  Gc,i  contains 
{mf/2)  X  (n//2)  cells  (full  coarsening).  Let  Tj  and  Q^j  denote  the  operators  that  transfer 
the  fine-grid  solution  and  residual  to  the  coarse  grid.  Since  any  coarse-grid  point  belongs 
also  to  the  fine  grid,  we  have  simple  injection  for  the  solution  transfer  operator.  Thus, 

T}'^Wf  =  Wf,  (3.20) 

To  ensure  the  conservation  property  for  the  residual  transfer,  we  defin^ 

(3.21) 
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where  fix  and  fiy  are  the  standard  averaging  operators  for  the  x  and  y  directions.  The 
symbol  of  this  operator  is  given  by 


=  (1  -f  cos^^)(l  -H  cos^,,). 

(3.22) 

If  we  now  apply  the  p-stage  scheme  ,  we  obtain  at  the  corresponding  point  of  the  Gc,i 
problem 

(»c,l)*'’  =(Wo,l)“  -a|<r.,l  [(ilc.l)''  ‘'  +  Fe,i  , 

=  (■Oci)*'". 

where  the  forcing  function 

1  =  1,2,-  -^p 

(3.23) 

PcA  =  Q‘,(Rf)*  -  {RcJ°\ 

and 

(R,)+  =  A{l,C,{w,)* 

By  successive  substitution  in  (3.23),  one  can  easily  show  that 

(3.24) 

Swc,i  =  -Fc,iSc,\(Tc,iQy^{Rfy 

or 

(3.25) 

^n)c,i  —  ~Fc,iSc,i(^c,iQ y 

with 

(3.26) 

(iD|)+  =  (w/)"  -1-  Swf. 

Since 

(3.27) 

Swf  =  -FZ{wf)^, 

then 

(3.28) 

(«)/)+  =  gfiwfY-, 

where  the  fine-grid  amplification  factor  is 

(3.29) 

1 

1— 1 

II 

(3.30) 

Moreover,  in  this  case,  where  we  are  considering  full  coarsening,  the  fine-grid  approximation 
is 


-f  6wf  +  6wc,u 


and  thus. 


(w,r+'  =  [l  - 


(3.31) 
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with  Qc  denoting  the  coarse-grid  amplification  factor. 

We  now  apply  the  multilevel  tmiform  analysis  to  two  semicoarsening  multigrid  strate¬ 
gies  (see  Figure  1).  Let  Gc,i  and  Gc,2  be  two  coarse  grids  containing  (m//2)  x  n/  and 
rrif  X  {nfl2)  cells,  respectively.  In  the  case  of  semicoarsening  with  simple  averaging,  we 
express  the  Fourier  transform  of  the  update  for  the  fine-grid  solution  as 

=  (u)/)"  (3.32) 


where  wi 
given  by 


UJ2  —  1/2.  Then  the  amplification  factor  associated  with  the  coarse  grids  is 


1 


ffc  —  1  ~  r  Fc,lSc,l<^c,lQ  f  ^  Zf -k- Fc, 2^0,2^  c,2Q  f  ^ 


c,2  -1 


(3.33) 


With  sequential  semicoarsening,  we  can  use  weightings  of  1.0  for  the  coarse-grid  corrections 
by  improving  the  estimate  of  the  initial  solution  on  either  grid  Gc,2»  as  indicated  in  Figure 
1,  or  Gc,i.  The  coarse-grid  correction  Sibc,i  is  given  by  (3.26),  and  the  other  one  Swc,2  is 
given  by 


SWc,2  =  -Ft 


c,2 


Sc,2<^c,2Q}'^  ~  ^c,2^c,lSc,\(Tc,lQy  ^ Z/jf/(uj/)".  (3.34) 


By  substituting  (3.26)  and  (3.34)  into  (3.32),  one  can  compute  the  amplification  factor 

Before  we  applied  the  two-dimensional  multigrid  analysis  just  discussed,  we  examined 
the  stability  properties  for  a  large  number  of  multistage  time-stepping  schemes,  including 
the  schemes  published  in  [19,20].  We  performed  this  preliminary  study  with  the  one¬ 
dimensional  advection  equation  for  a  single  grid.  Figures  2  and  3  show  two  examples  from 
this  study.  The  dash  line  represents  the  locus  of  the  symbol  of  the  difference  operator, 
and  it  must  lie  inside  the  absolute  stability  curve.  The  dissipative  term  for  first-order 
upwind  differencing  is  evaluated  three  times  for  both  the  three-  and  five-stage  schemes. 
In  the  case  of  the  (5,3)  scheme,  the  dissipation  is  computed  on  the  first,  third,  and  fifth 
stages,  and  the  weights  are  those  given  in  (3.14).  While  the  (3,3)  scheme  exhibits  fairly 
good  high-frequency  damping,  there  is  a  substantial  improvement  in  damping  with  the 
(5,3)  scheme,  at  the  expense  of  a  little  extra  work.  The  coefficients  for  the  (5,3)  scheme 
are  those  determined  by  Tal  [21]. 

In  Figures  4  and  5,  we  present  contours  of  the  amplification  factor  for  the 

following  cases;,  1)  single  grid,  2)  full  coarsening,  3)  semicoarjening  with  simple  averaging, 
4)  sequential  semicoarsening.  For  each  case,  cell  aspect  ratios  of  one  and  ten  are  shown. 
The  improvement  in  damping  with  the  two-level  schemes  is  evident.  There  is  greater 
compression  of  to  the  origin  with  the  semicoarsening  strategies  The  aspect  ratio 

(A)  of  ten  contours  indicate  that  the  one-dimensional  behavior  of  the  driving  scheme  in 
the  T]  direction  is  recovered  for  large  A.  They  also  show  that  modes  associated  with  the  | 
diix  cion  are  damped  much  better  with  the  semicoarsening  schemes. 

According  to  the  current  analysis,  where  the  nodal  point  of  interest  is  assumed  to  be 
common  to  all  meshes  being  considered,  the  semicoarsening  scheme  with  selective  averaging 
has  the  same  damping  behavior  as  the  one  with  simple  averaging.  In  practice,  however, 
we  observe  the  expected  improvement  in  damping  when  using  selective  rather  than  simple 
averaging.  This  will  be  demonstrated  later. 
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4  Governing  Flow  Equations 

Let  p,u,v,p,T,E^  and  H  denote,  respectively,  the  nondimensional  values  of  density, 
velocity  components  in  the  x  and  y  Cartesian  directions,  pressure,  temperature,  specific 
total  internal  energy,  and  specific  total  enthalpy.  In  addition,  let  Ci  and  be  unit  vectors 
of  the  Cartesian  coordinate  system  (x,  y),  and  let  n  be  a  unit  vector  normal  to  the  surface 

5  enclosing  a  volume  V.  Then  the  two-dimensional,  unsteady  Navier-Stokes  equations, 
neglecting  body  forces  and  heat  sources,  can  be  written  in  conservative  form  in  a  Cartesian 
coordinate  system  as 

^  IJv 

where  the  solution  vector  W  and  the  tensors  Tc,  Tv  are  defined  as 


■  p  ■ 

pUCx  -t-  pVBy 

^  0 

pu 

{pii'^  +  p)ex  -f-  puvey 

11 

1 

Ox^x  "t"  Txy^y 

w  = 

pv 

II 

puvGx  +  (pv^  +  p)ey 

Txy^x  +  O^^y 

pE 

puHex  pvHey 

m 

{uOx  +  VTxy  — 

+  {xiTxy+Vay-qy)ey 

with 


Ox  — 
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.  ( du  dv\ 


H=^E+^, 

P 


The  scaling  factor  Re  ^  =  ,yjMRe~^,  with  M  and  Re  representing  the  Mach  and 
Reynolds  numbers,  respectively.  In  this  paper,  the  working  fluid  is  air,  and  it  is  assumed 
to  be  thermally  and  calorically  perfect.  That  is,  the  equation  of  state  is 


/.  =  (7~l)/.(£-(uHt>^)/2),  T=p/p  (4.2) 
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The  quantities  fi  and  A  are  the  first  c..id  second  coefficients  of  viscosity,  respectively,  and 
A  is  taken  to  be  (Stokes  hypothesis).  Sutherland’s  law  is  used  to  determine  the 
molecular  viscosity  coefficient  fi.  The  coefficient  of  thermal  conductivity  (k)  is  evaluated 
using  the  constant  Prandtl  number  assumption.  The  effect  of  turbulence  is  taken  into 
account  by  using  the  eddy-viscosity  hypothesis.  In  the  present  work,  the  turbulence  model 
of  Baldwin  and  Lomax  [22]  is  used. 


5  Spatial  Discretization 


The  numerical  approximation  of  (4.1)  follows  the  method  of  lines,  which  decouples  the 
discretization  in  space  and  time.  The  physical  domain  around  the  aerodynamic  body  is 
divided  into  quadrilateral  cells  by  the  generation  of  a  body-fitted  grid.  The  discrete  values 
of  the  flow  quantities  are  located  at  the  vertices  of  the  mesh  cells.  For  the  flux  calculation, 
an  auxiliary  grid  is  used,  which  is  defined  by  connecting  the  cell  centers  of  the  original  cells 
(see  Figure  6).  The  integral  equation  (4.1)  is  approximated  by  the  spatial  discretization, 
yielding 

d  1 


(5.1) 


where  Vij  denotes  the  area  of  the  control  volume  surrounding  the  grid  node  (i,  j).  (Qc)ij 
and  (Qti)ij  represent  approximations  of  the  convective  flux  and  viscous  flux,  respectively. 
The  viscous  fluxes  are  approximated  by  central  differences  using  a  local  transformation 
from  Cartesian  coordinates  to  the  curvilinear  coordinates  [4],  In  the  following  some  de¬ 
tails  of  the  upwind  discretization  of  the  convective  terms  are  discussed.  The  inviscid  flux 
through  the  interface,  (i  4- 1/2,  j),  is  evaluated  as 


(Qc), 4.1/2, j  -  2  [('^c)i,j  +  (<?>)i-H,>]  -  S, 4.1/2, j  -f-  -Ri+ii2,j^,+  i/2,y  (5.2) 

Here,  S, 4.1/2, j  is  the  surface  vector  of  face  {i  +  ll2,j),  and  R  is  the  right  eigenvector 
matrix  of  the  flux  Jacobian  in  transformed  space.  Equation  (5.2)  separates  the  inviscid 
numerical  flux  into  the  sum  of  an  averaged  term  corresponding  to  central  differencing  and 
a  dissipative  term,  which  adapts  the  discretization  stencil  in  accordance  with  local  wave 
pr  '■nation.  The  flux  function  ^  is  based  on  the  second-order  accurate  upwind  TVD 
schen...  of  Yee  and  Harten  [23].  Here,  we  do  not  repeat  the  formulation  for  $,  but  simply 
indicate  some  important  details  in  the  present  numerical  evaluation  of  Second-order 
accuracy  is  obtained  with  the  limiter 


Q^i-l/2,j(Q|+l/2,j^  +  e)  +  Q, 4-1/2, jjQt-l/Z.j^  4-  g) 

<^i-i,/2,j^  +  +  2e 


(5.3) 


where  a  is  the  first  difference  of  the  characteristic  variable,  and  e  is  a  small 

constant  to  prevent  division  by  zero.  The  flow  quantities  at  the  face  {i  +  1/2, j)  are 
evaluated  by  Roe’s  averaging  procedure  [24].  Note  that  the  scheme  is  identical  to  Roe’s 
first-order  flux  difference  splitting  for  g  —  0,  Since  Roe’s  scheme  may  violate  the  entropy 
condition  when  the  eigenvalues  of  the  flux  Jacobians  vanish,  the  eigenvalues  are  modified 


11 


using  an  entropy  function  This  entropy  correction  has  to  be  carefully  designed  for 
viscous  flow  calculations.  The  shear  layers  along  solid  walls  are  numerically  smeared,  if 
an  entropy  correction  is  applied  to  the  eigenvalues  associated  with  the  convective  waves. 
On  the  other  hand,  if  cells  with  high  aspect  ratios  are  present,  additional  support  for 
damping  in  the  direction  of  the  long  side  of  a  cell  is  needed  in  regions  of  low  velocities, 
such  as  stagnation  points.  Therefore,  the  correction  of  the  eigenvalues  ’J*  is  constructed  as 
a  function  of  the  cell  aspect  ratio.  For  the  convective  waves,  the  function 

+(A  otherwise, 


and  for  the  acoustic  waves 


$(A(')  = 


|A<'|. 


if  |A{'|  > 
otherwise. 


(5,5) 


In  equations  (5.4)  and  (5.5),  represents  the  /-th  eigenvalue  of  the  transformed  flux 
Jacobian  in  the  ^-direction.  The  entropy  parameter,  S,  is  given  according  to  Muller  [25] 


^  =  a«(i  +  (A,/A^)"), 


(5.6) 


where  A^,  A,,  are  the  spectral  radii  of  the  flux  Jacobians  in  the  ^  and  17  directions,  0  <  w  <  1, 
and  0.1  <  ^  <  0.5.  If  U  =  ue*  +  ucy,  and  are  directed  areas  associated  with  the  ^ 
and  T)  directions,  and  c  is  the  speed  of  sound,  then 

A|  =  |U  ■  S«|  +  cIS^I,  A,  =  |U  .  S,|  +  c|S,|.  (5.7) 

The  blending  coefficient,  /?,  accounts  for  the  cell  aspect  ratio.  It  is  given  as 


=  max(l  -  A^/A,,,/«).  (5.8) 

It  is  shown  below  that  k  should  be  zero  for  accurate  computations  of  shear  layers.  Fur¬ 
thermore,  we  will  demonstrate  that  a  wide  range  of  flow  problems  can  be  solved  accurately 
with  a  single  set  of  parameters,  6  —  0.25,  u  =  0.3,  and  «  =  0. 

6  Time-Stepping  Scheme  for  Hypersonic  Flows 

The  basic  elements  of  the  time-stepping  scheme  have  been  outlined  in  Section  3, 
and  they  are  not  repeated  here.  In  the  following  sections  emphasis  is  placed  on  recent 
improvements  to  enhance  robustness  for  hypersonic  flow  calculations. 

6.1  Multistage  Scheme  for  the  Fine  and  Coarse  Meshes 

Consistent  with  the  results  for  the  advection  equation  in  Section  3,  we  have  observed 
the  need  to  pair  spatial  discretization  and  particular  time-stepping  schemes  for  the  solution 
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of  the  Navier-Stokes  equation.  The  most  robust  choice  of  spatial  discretizations  found  to 
this  point  is  to  use  the  second-order  upwind  scheme  of  Section  5  on  the  fine  meshes  and  to 
set  the  limiter  to  zero  everywhere  on  the  coarse  meshes.  An  alternative  choice  taken  in  [6,8] 
is  to  use  scalar  second-difference  dissipation  terms  on  the  coarse  meshes.  This  turned  out  to 
be  less  robust  because  the  second  differences  are  less  diffusive  with  respect  to  the  acoustic 
modes;  also,  the  central-difference  scheme  allows  waves  to  travel  upstream  in  supersonic 
flow.  As  indicated  previously,  a  five-stage  scheme  with  three  evaluations  of  dissipation  and 
the  coefficients  of  Tai  [21]  is  used  for  time  advancement.  The  large  viscous  stability  limit 
of  this  scheme  is  shown  in  Figure  3.  Disturbances  are  most  effectively  expelled  out  of  the 
computational  domain  by  using  local  time  stepping  and  implicit  residual  smoothing  [5,6].. 
The  smoothing  of  the  residuals  allows  a  Courant  number  ratio  (CFLfCFL*)  as  high  as 
2.5  when  CFL*  =  2.3,  which  is  roughly  the  stability  limit  of  the  explicit  scheme.  The 
time  step  is  determined  by  the  spectral  radii  of  the  inviscid  flux  Jacobians  in  the  different 
coordinate  directions,  and  A,,,  as 

A(  =  (6.1) 

-t-  A,, 

In  order  to  stabilize  the  schemes  in  regions  where  the  viscous  stability  limit  is  more 
restrictive  than  the  inviscid  limit,  the  coefficients  of  the  implicit  residual  smoothing  oper¬ 
ator  are  locally  increased,  as  outlined  in  [6,8].  At  strong  shocks,  however,  high  Courant 
numbers  are  not  appropriate.  Consequently,  an  adaptive  time  step  is  employed.  By  using 
the  nondimensional  second  difference  of  the  pressure  as  a  switch,  the  value  of  CFL  is 
locally  reduced  to  about  two  at  the  shock. 

6.2  Multigrid  Schemes 

For  the  numerical  solution  of  the  Navier-Stokes  equations,  the  two-level  strategies 
presented  in  Figure  1  are  extended  to  multilevel  schemes,  as  displayed  in  Figure  7.,  The 
only  differences  between  the  two-level  schemes  and  the  multilevel  schemes  occur  in  the 
restriction  process.  Whenever  two  “down”  arrows  meet  at  a  coarse  mesh,  averaging  is 
used  to  obtain  the  restricted  variable.  The  multilevel  arrangement  of  the  coarse  meshes, 
shown  in  Figure  7  b),  was  first  given  by  Mulder  [17],  who  used  semicoarsening  in  order 
to  solve  the  flow  alignment  problem.  Suitable  coordinate  meshes  for  thin  boundary  layers 
exhibit  mostly  cells  with  high  aspect  ratios  in  the  surface-aligned  direction.  Figure  8 
displays  further  variants  of  semicoarsening  for  these  situations  which  are  computationally 
cheaper  than  the  semicoarsening  schemes  shown  in  Figure  7. 

One  may  notice  that  the  central  restriction  and  prolongation  operators  discussed  in 
Section  2  allow  for  upstream  propagation  of  disturbances  in  supersonic  flow.;  Further¬ 
more,  the  corrections  given  by  the  standard  multigrid  scheme  near  strong  shocks  lead  to 
divergence  of  the  calculation,  particularly  during  the  initial  part  of  the  transient  phase. 
Therefore,  the  restriction  operator  is  damped  by  using 

i?,,j  =  max(l  -  (6.2) 

where  Ri^j  is  the  standard  restriction  operator,  and  is  a  switch  to  detect  strong 

shocks, 

=  ^^"Vax(i/,,i/j+i,i/,_i,z/j,i/j+i,r/j_i),  (6.3) 
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= 


(6.4) 


Pi-1,3  d*  Pi-1-1, j 

Pt,3-1  -  ^Pi,j  +  Pi,3+l 

Pi-1,3  d"  2pi,j  +  Pi+l,3 

,  Uj  — 

Pi,3-l+^Pi,3  +  Pi, J+1 

The  damping  coefficient,  k^"'\  is  given  a  value  of  about  one  in  the  start-up  phase  of  the 
multigrid  process  and  is  decreased  to  a  value  of  about  0.4  at  later  cycle  numbers  in  order  to 
allow  for  good  asymptotic  convergence  rates.;  This  is  in  line  with  the  restriction  damping  of 
Koren  and  Hemker  [14],  who  based  their  damping  coefficients  on  a  more  physical  analysis. 

A  fixed  V-type  cycle  with  time  stepping  on  the  way  down  is  used  to  execute  the  multi¬ 
grid  strategies  described  above.;  The  robustness  of  the  overall  scheme  is  much  improved 
by  smoothing  the  resultant  coarse-mesh  corrections  before  they  are  passed  to  the  finest 
mesh.  The  smoothing  reduces  high-frequency  oscillations  introduced  by  the  linear  interpo¬ 
lation  of  the  coarse-mesh  corrections.  The  factored  scheme  equation  (3.16)  with  constant 
coefficients  around  0.1  is  used  for  this  smoothing.  Also,  the  application  of  Full  Multigrid 
(FMG)  provides  a  well-conditioned  starting  solution  for  the  finest  mesh  being  considered. 

7  Numerical  Results 

Three  different  hypersonic  flow  cases  axe  used  to  assess  the  capabilities  of  the  multigrid 
schemes.  These  are  laminar  Mach  10  (M  =  10)  flow  over  a  compression  ramp,  turbulent 
flow  over  a  slender  forebody  at  high  Reynolds  numbers,  and  laminar  flow  over  an  airfoil  at 
high  Mach  number  and  high  angle  of  attack.  Table  1  gives  a  summary  of  the  geometries 
and  the  flow  parameters  of  the  test  cases.  In  this  table,  Tjn/  is  the  dimensional  free- 
stream  temperature,  and  Tu,  is  the  specified  wall  temperature.  Also,  the  finest  grid  used 
for  each  flow  computation  is  characterized  by  the  streamwise  and  normal  leading-edge 
spacings  (As/e,  An/e),  along  with  the  normal  spacing  (Au/e)  at  the  end  or  trailing  edge  of 
a  geometry. 

The  flow  over  the  compression  ramp  is  identical  to  Case  3.2  of  the  Workshop  on 
Hypersonic  Flows  for  Reentry  Problems,  Part  II,  held  in  Antibes,  France,  1991.  This 
allows  comparisons  with  the  performance  of  other  computational  methods  published  in 
[26] .  Figure  9  displays  the  coordinate  mesh  generated  for  this  test  case.  The  low  Reynolds 
number  allows  for  a  mesh  with  moderate  aspect  ratios  between  5  and  50  near  the  wall. 
The  129  X  81  mesh  is  successively  coarsened  down  to  9  x  6,  which  yields  9  grid  levels  with 
semicoarsening  and  5  levels  with  full  coarsening.  It  is  expected  that  the  semicoarsening 
strategy  should  eliminate  most  of  the  stiffness  associated  with  aspect  ratio  .  The  converged 
flow  solution  is  shown  in  Figure  10,  The  computed  extent  of  separation  in  the  corner  is 
somewhat  smaller  for  the  coarse  mesh  than  for  the  fine  mesh.  Note  that  the  result  of  the 
fine  mesh  agrees  very  well  with  grid- converged  computations  published  in  [13], 

In  the  next  figures,  we  analyze  the  performance  of  the  different  multigrid  schemes. 
For  this  purpose,  computations  are  started  from  a  solution  which  was  converged  to  about 
plotting  accuracy.  Figure  11  compares  the  different  schemes  of  Figure  7,  The  numbers 
indicate  the  final  convergence  rate  (r)  of  the  schemes  and  the  rate  of  data  processing 
{RDP)  on  a  CRAY-YMP  to  advance  one  grid  point  by  one  multigrid  cycle.  It  is  seen  that 
the  sequential  semicoarsening  scheme  (Figure  7  c))  gives  by  far  the  best  convergence  rate. 
For  this  scheme,  the  effect  of  the  modifications  shown  in  Figure  8  is  investigated  in  Figure 
12.  One  finds  that  the  meshes  obtained  by  full  coarsening  and  by  semicoarsening  in  the 
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direction  normal  to  the  wall  are  both  important  to  achieve  good  convergence  rates.  From 
Figures  11-12,  we  conclude  that  semicoarsening  with  a  selected  number  of  coarse  meshes 
is  most  effective  for  this  flow  problem;  however,  full  coarsening  does  a  surprisingly  good 
job  because  of  its  low  work  count.  Figures  13-14  show  improvements  which  may  be  gained 
by  using  more  than  a  single  time  step  on  the  coarse  meshes.  The  full  coarsening  scheme. 
Figure  13,  gives  only  marginal  gains  when  using  more  than  two  time  steps  on  the  coarse 
meshes.  The  sequential  semicoarsening  scheme.  Figure  14,  gives  an  initially  improved 
rate,  whereas  the  final  rate  is  not  affected  by  more  work  done  on  the  coarse  meshes.  It  is 
thought  that  the  capabilities  of  the  multigrid  approach  are  put  to  full  use  for  this  test  case. 
Further  improvements  are  forseen  only  if  the  remaining  stiffness  in  the  discrete  equations, 
that  is  the  differences  in  the  characteristic  speeds  of  acoustic  and  convective  waves,  can 
be  overcome  by  some  proper  means.  A  comparison  between  single  mesh  and  multigrid 
computations  is  given  in  Figure  15.  We  find  that  the  multigrid  scheme  with  sequential 
semicoarsening  converges  within  one  tenth  of  the  computing  time  required  for  the  single 
mesh  scheme.; 

The  grid  generated  for  the  ramp  flow  is  well  suited  to  study  the  grid-alignment  problem 
which  occurs  for  inviscid  flow  over  the  ranip.  Figure  16  shows  convergence  histories  of 
various  schemes  obtained  by  using  a  slip- wall  boundary  condition  and  omitting  the  viscous 
terms  in  the  governing  equations.  Generally,  convergence  is  worse  for  the  inviscid  flow 
because  the  convective  eigenvalues  of  the  flux  Jacobian  in  the  normal  direction  are  zero 
for  most  of  the  grid  points.  The  second-order  solution  does  not  converge,  regardless  which 
multigrid  strategy  we  use..  With  the  flux  limiter,  equation  (5.3),  set  to  zero  everywhere  on 
the  fine  mesh,  the  solution  converges,  provided  we  use  all  the  coarse  meshes  introduced  by 
the  semicoarsening  approach..  Note  that  this  result  is  identical  to  the  findings  of  Mulder 
[27].  We  have  to  conclude  that  the  present  multigrid  schemes  cannot  cope  with  the  grid- 
alignment  problem  in  the  context  of  high-order  spatial  discretizations. 

The  flow  over  a  slender  forebody  is  chosen  to  represent  a  generic  configuration  corre¬ 
sponding  to  a  high-speed  civil  transport  aircraft  or  an  air-breathing  space  transportation 
system  with  low  wave  drag.  The  high  Reynolds  numbers  yield  thin  boundary  layers,  which 
can  only  be  resolved  with  highly  clustered  coordinate  meshes  and  large  aspect  ratio  cells. 
The  mesh  used  for  the  present  investigations  is  di.splayed  in  Figure  17.  The  cells  near  the 
wall  have  aspect  ratios  up  to  25000.  The  flow  computations  were  done  with  fixed  transi¬ 
tion  at  2  percent  chord  and  with  the  assumption  of  an  adiabatic  wall..  Figures  18-20  show 
the  solution  obtained  on  three  successively  refined  meshes..  Both  the  distributions  of  skin 
friction  and  wall  temperature  are  accurately  computed,  even  with  just  25  points  in  the 
normal  direction.  The  effect  of  numerical  dissipation  introduced  by  the  entropy  correction 
function  4*,  equations  (5.4)-(5.8),  is  shown  in  Figure  21.  If  the  value  of  a;  is  increased  to 
2/3,  which  is  a  typical  value  used  in  central-difference  codes  [2,4],  the  wall  temperature  is 
badly  reproduced  on  the  coarse  mesh.  Furthermore,  if  we  introduce  numerical  dissipation 
for  the  convective  waves  in  the  direction  of  the  short  side  of  the  cells  by  setting  k  >  0, 
the  shear  layers  are  numerically  smeared;  and  hence,  the  wall  temperature  is  adversely 
affected..  We  notice  that  a  value  of  k  =  0.1,  which  again  is  a  representative  dissipation 
level  of  current  central- difference  codes,  yields  reasonable  solutions.  However,  the  conver¬ 
gence  behavior  is  not  changed  much  with  this  relatively  low  level  of  numerical  dissipation. 
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Hence,  all  the  results  presented  for  the  three  test  cases  were  obtained  with  a  single  set  of 
parameters,  6  =  0.25,  u  =  0.3,  and  /c  =  0. 

Next  we  investigate  the  convergence  behavior  of  the  multigrid  schemes.,  The  fine  mesh 
with  257  X  97  points  allows  11  grid  levels  to  be  used  with  semicoarsening.  First  we  notice 
that  it  is  not  possible  to  run  the  full  diamond-shape  tree  of  coarse  meshes.  Obviously,  the 
time-stepping  scheme  is  not  well  suited  to  handle  the  extreme  aspect  ratios  which  occur 
on  the  17  x  97  mesh,  for  example.  Using  the  proper  half  of  the  diamond,  which  includes 
the  meshes  with  relatively  low  aspect  ratios,  the  numerical  solution  converges.  Figure  22 
displays  a  comparison  of  the  different  multigrid  strategies.  The  computations  are  started 
from  a  preconverged  solution..  The  computer  time  to  update  a  grid  point  by  one  multigrid 
cycle,  RDP,  and  the  final  convergence  rate,  r,  are  also  included.  Again,  the  scheme  with 
sequential  semicoarsening  converges  best.  The  differences  between  the  multigrid  schemes 
for  this  case,  having  cells  with  very  high  aspect  ratios,  are  larger  than  for  the  ramp  flow. 
The  final  convergence  rate  of  the  scheme  with  sequential  semicoarsening  is  15  times  better 
than  the  rate  with  full  coarsening.  A  comparison  of  the  performance  for  the  complete 
FMG  process  is  given  in  Figure  23.  The  sequential  semicoarsening  scheme  takes  194  cycles 
and  570s  CPU  on  a  CRAY-YMP  to  reduce  the  averaged  residuals  to  10“^  on  the  fine 
mesh.  The  scheme  with  full  coarsening  takes  1024  cycles  and  1430s,  and  the  single  mesh 
code  takes  7762  time  steps  and  6190s  to  achieve  the  same  convergence  level.  Note  that 
residuals  of  10“^  correspond  to  a  solut  on  which  is  converged  within  plotting  accuracy.  If 
we  compared  computer  times  to  reach  >ower  levels  of  residuals  instead,  the  results  would 
have  been  even  better  for  the  multigrid  scheme  with  semicoarsening. 

Laminar  flow  over  an  airfoil  at  M  =  25  and  a  =  30  degrees  is  chosen  as  a  final 
test  case  in  order  to  demonstrate  that  the  multigrid  method  used  here  can  handle  very 
strong  shock  waves  and  highly  expanded  flow.  Figure  24  shows  the  257  x  81  mesh  .  The 
numerical  flow  solution,  which  is  represented  in  Figures  25-28,  features  a  large  separated 
flow  region  with  two  distinct  vortices.  The  resulting  shear  layers  are  not  well  aligned 
with  the  coordinate  mesh;  and  hence,  considerable  numerical  smearing  is  expected  in 
those  regions.  The  difficulties  in  resolving  this  highly  separated  flow  are  illustrated  by  a 
comparison  of  the  results  obtained  from  meshes  with  different  grid  densities.  Obviously, 
the  mesh  with  129  x  41  points  is  still  too  coarse  to  resolve  the  separated  flow  region.  The 
convergence  history  of  the  sequential  semicoarsening  scheme  is  shown  in  Figure  29.  The 
residual  drops  8  orders  of  magnitude  within  400  multigrid  cycles.  Note  that  the  largest 
values  of  dp/dt  are  located  at  the  shock  wave.  The  solution  converges  to  plotting  accuracy 
within  100  multigrid  cycles. 

8  Conclusions 

New  multigrid  schemes  for  hypersonic  flow  computations  have  been  investigated.  The 
basic  solution  algorithm  employs  upwind  discretization  and  explicit  multistage  time  step¬ 
ping.  Various  multigrid  schemes  with  semicoarsening  are  introduced  in  order  to  overcome 
stiffness  resulting  from  the  mesh  cells  with  high  aspect  ratio  which  are  necessary  to  resolve 
viscous  flows.  The  basic  components  of  the  algorithm  are  examined  with  Fourier  stabiliy 
analysis  applied  to  the  two-dimensional  advection  equation.  Both  the  results  of  the  Fourier 
analysis  and  the  computations  of  high  Reynolds  number  flows  suggest  that  the  semicoars- 
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ening  approach  is  effective.  For  the  first  time,  convergence  rates  for  hypersonic  viscous 
flows  axe  shown  which  are  similar  or  even  better  than  those  previously  published  for  the 
transonic  regime  [3,4].,  Further  work  is  required  to  make  the  computational  scheme  less 
expensive.  This  is  especially  true  for  the  coarse  meshes  used  within  the  semicoarsening 
approach,  which  make  up  the  major  portion  of  the  overall  work  count  of  the  scheme.  Fur¬ 
ther  improvements  of  convergence  rates  seem  possible  if  stiffness  arising  from  the  difference 
of  characteristic  speeds  of  acoustic  and  convective  waves  can  be  overcome.  For  this  pur¬ 
pose,  new  techniques  such  as  preconditioning  of  the  flow  equations  [28,29]  or  characteristic 
residual  smoothing  [12]  seem  to  hold  promise. 
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Table  1:  Flow  parameters  and  geometric  parameters  for  the  test  cases 


(in,n) 

'  1 
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(in/2,n/2) 
a)  Full  coarsening 

(ro,n) 


c)  Sequential  semicoarsening 

Figure  1: 


(m,n) 


b)  Semicoarsening  with  simple  averaging 


(m,  n) 


Semicoarsening  with  selective  averaging 


d) 


Two-level  multigrid  schemes  investigated  in  the  present  work 
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(a)  Stability  curves  with  implicit  residual  smoothing  (CFL  =  2.6,  /?  =  0.75) 
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(b)  Amplification  factor  with  implicit  residual  smoothing  (CFL  =  2.6,  ^  =  0.75) 

Figure  2  Stability  plots  for  3-stage  Runge-Kutta  scheme  with 
first-order  upwind  approximation  (coefficients :  0.105,  0.325,  1.0) 
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(a)  Stability  curves  with  implicit  residual  smoothing  (C?L  =  5.0,  ^  -  1.0) 
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(b)  Amplification  factor  with  implicit  residual  smoothing  (CFL  =  5.0,  /?  =  1.0) 

Figure  3  Stability  plots  for  5-stage  Runge-Kutta  scheme  with  first-order  upwind  approximation 
and  3  evaluations  of  dissipation  (coefficients ;  0.2742,  0.2067,  0.6020,  0.5142,  1.0) 
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Figure  4  Contour  plots  of  amplification  factor  for  5-stage  Runge-Kutta 
scheme  with  first-order  upwind  approximation  and  3  evaluations 
of  dissipation  (coefficients  0.2742,  0.2067,  0.5020,  0.5142,  1.0) 


(c)  Two  levels,  semicoarsening  with  simple  averaging, 
weights  =  0.5,  A  =  1  (CFL  =  5.0,  CFL*  =  2.4) 


Figure  4  Contour  plots  of  amplification  factor  for  5-stage  Runge-Kutta 
scheme  with  first-order  upwind  approximation  and  3  evaluations 
of  dissipation  (coefficients  :  0.2742,  0.2067,  0.5020,  0.5142,  1.0) 


(b)  Two  levels,  full  coarsening,  A  =  10  (CFL  =  5.0,  CFL  =  2.4) 


Figure  5  Contour  plots  of  amplification  factor  for  5-stage  Runge-Kutta 
scheme  with  first-order  upwind  approximation  and  3  evaluations 
of  dissipation  (coefficients  ;  0.2742,  0.2067,  0.5020,  0.5142,  1.0) 


(c)  Two  levels,  semicoarsening  with  simple  averaging, 
weights  =  0.5,  A  =  10  (CFL  =  5.0,  CFL*  =  2.4) 


(d)  Two  levels,  sequential  semicoarsening,  weights  =  1.0,  A  =  10  (CFL  =  5.0,  CFL*  =  2.4) 


Figure  5  Contour  plots  of  amplification  factor  for  5 -stage  Runge-Kutta 
scheme  with  first-order  upwind  approximation  and  3  evaluations 
of  dissipation  (coefficients  :  0.2742,  0.2067,  0.5020,  0.5142,  1.0) 
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Figure  6:  Control  volume  for  nodal-point  scheme 


c)  Sequential  semi  coarsening  d)  Semicoarsening  with  selective  averaging 


Figure  7:  Multilevel  schemes 
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Figure  8: 


Semicoarsening  with  selection  of  coarse  meshes 
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b)  Pressure  coefiicient  c)  Skin  friction 

Figure  10:  Flow  solution  for  ramp-flow  problem 
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Figure  11:  Influence  of  multigrid  strategies  on  convergence  for  ramp-flow  problem 


Figure  12:  Influence  of  selected  coarse  meshes  on  convergence  for  ramp-flow  problem 
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Figure  13:  Influence  of  multiple  time  stepping  on  coarse  meshes  for  multigrid  with  full 
coarsening 


Figure  14:  Influence  of  multiple  time  stepping  on  coarse  meshes  for  multigrid  with  sequen¬ 
tial  semicoarsening 
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Figure  15:  Convergence  of  ramp  flow  with  single-mesh  time  stepping  and  multigrid  with 
sequential  semicoarsening 


Figure  16:  Convergence  histories  for  flow-alignment  problem 
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Figure  19:  Distribution  of  skin  friction  along  forebody 


Figure  20:  Distribution  of  adiabatic  wall  temperature  along  forebody 


b)  Influence  of  the  cut-off  value,  k,  applied  to  the  convective  waves 


Figure  21:  Influence  of  the  entropy  correction  coefficients  on  the  adiabatic  wall  temperature 
along  forebody 


multigrid  cycles 

Figure  22:  Influence  of  multigrid  strategy  on  convergence  for  forebody,  mesh  256x96 
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Figure  23:  Convergence  histories  for  single-mesh  time  stepping  and  multigrid  with  sequen¬ 
tial  semicoarsening 
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Figure  24:  Coordinate  mesh  for  NACA  0012  Figure  25:  Mach  contours  for  NACA  0012 
with  256x80  cells  airfoil  at  M  =  25,  a  =  30® 


Figure  26:  Streamlines  in  separated-flow  region 
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Figure  27:  Distribution  of  skin  friction  along  airfoil 


Figure  28:  Distribution  of  Stanton  number  along  airfoil 
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Figure  29:  Convergence  history  for  flow  around  NACA  0012  airfoil  at  M  =  25,  a  =  30®. 
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